Initialize

   if(any(grepl("package:dvutils", search()))) detach("package:dvutils", unload=TRUE)
   library(dvutils)
   library(stringr)
   library(ggplot2)
   library(dplyr)
   library(tidyr)
   library(tibble)
   library(readr)
   library(maps)
   library(glue)
   
   library(phyloseq)

   # library(pr2database)
   # data("pr2")

## Warning: package 'knitr' was built under R version 3.5.3

Get the Antarctica ASVs

• Only use asv for which supergroup_boot >= 90

rm(list = ls())
sample_type <- c("18S filter", "16S plastid", "18S sort")

asv <- list()
ps <- list()
long <- list()
# Export as specific data set as a phyloseq file

i <- 16
for (one_sample_type in sample_type) {
    asv[[one_sample_type]] <- metapr2_export_asv(dataset_id_selected = i, export_phyloseq = FALSE, 
        export_xls = FALSE, boot_min = 90, boot_level = supergroup_boot, directory = "../dada2/")
    
    ps[[one_sample_type]] <- asv[[one_sample_type]][["ps"]]
    # Label samples with date
    sample_data(ps[[one_sample_type]])$sample_label <- as.character(sample_data(ps[[one_sample_type]])$date)
    
    # Add TFF to sample name sample_data(ps[[one_sample_type]])$sample_label <-
    # str_c(sample_data(ps[[one_sample_type]])$metadata_code,
    # str_replace_na(sample_data(ps[[one_sample_type]])$sample_concentration,
    # ''), sep = '.' )
    
    print(glue("Phyloseq - {one_sample_type}"))
    print(ps[[one_sample_type]])
    cat("============================\n")
    
    i <- i + 1
    
}

Phyloseq - 18S filter
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 2310 taxa and 123 samples ]
sample_data() Sample Data:       [ 123 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 2310 taxa by 8 taxonomic ranks ]
============================
Phyloseq - 16S plastid
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 442 taxa and 103 samples ]
sample_data() Sample Data:       [ 103 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 442 taxa by 8 taxonomic ranks ]
============================
Phyloseq - 18S sort
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 385 taxa and 60 samples ]
sample_data() Sample Data:       [ 60 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 385 taxa by 8 taxonomic ranks ]
============================

sample_variables(ps[["18S filter"]])

 [1] "sample_id"                      "file_name"                     
 [3] "sample_name"                    "sample_code"                   
 [5] "metadata_code"                  "replicate"                     
 [7] "DNA_RNA"                        "fraction_name"                 
 [9] "fraction_min"                   "fraction_max"                  
[11] "sample_concentration"           "sample_sorted"                 
[13] "reads_total"                    "sample_remark"                 
[15] "metadata_id"                    "metadata_code_original"        
[17] "project"                        "cruise"                        
[19] "station_id"                     "station_id_num"                
[21] "year"                           "date"                          
[23] "time"                           "season"                        
[25] "depth_level"                    "depth"                         
[27] "substrate"                      "substrate_description"         
[29] "substrate_description_detailed" "experiment_name"               
[31] "experiment_time"                "experiment_time_unit"          
[33] "experiment_bottle"              "experiment_condition"          
[35] "latitude"                       "longitude"                     
[37] "site_name"                      "country"                       
[39] "oceanic_region"                 "bottom_depth"                  
[41] "temperature"                    "salinity"                      
[43] "pH"                             "O2"                            
[45] "fluorescence"                   "ice_coverage"                  
[47] "Chla"                           "NO2"                           
[49] "NO3"                            "PO4"                           
[51] "Si"                             "Chla_0.2_3 um"                 
[53] "bact_ml"                        "syn_ml"                        
[55] "peuk_ml"                        "neuk_ml"                       
[57] "crypto_ml"                      "virus_small_ml"                
[59] "virus_large_ml"                 "metadata_remark"               
[61] "sample_label"                  

rm(asv)

Filtration

• Only use Station 6 (not 14)
• Only surface considered (5 m)
• Do not consider TFF samples
• Only keep photosynthetic groups abd exclude dinoflagellates
• Very important, must remove taxa that are not present in the filtered samples

for (one_sample_type in sample_type) {
    
    ps[[one_sample_type]] <- ps[[one_sample_type]] %>% subset_taxa((division %in% 
        c("Chlorophyta", "Cryptophyta", "Rhodophyta", "Haptophyta", "Ochrophyta")) | 
        (class == "Filosa-Chlorarachnea")) %>% subset_samples(station_id == 
        "6") %>% subset_samples(depth_level == "surface") %>% subset_samples(is.na(sample_concentration)) %>% 
        filter_taxa(function(x) sum(x) > 0, TRUE)
    
    
    print(glue("Phyloseq - {one_sample_type}"))
    print(ps[[one_sample_type]])
    cat("============================\n")
    
}

Phyloseq - 18S filter
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 394 taxa and 50 samples ]
sample_data() Sample Data:       [ 50 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 394 taxa by 8 taxonomic ranks ]
============================
Phyloseq - 16S plastid
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 236 taxa and 40 samples ]
sample_data() Sample Data:       [ 40 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 236 taxa by 8 taxonomic ranks ]
============================
Phyloseq - 18S sort
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 114 taxa and 16 samples ]
sample_data() Sample Data:       [ 16 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 114 taxa by 8 taxonomic ranks ]
============================

Get the different fraction separately

fraction_filter <- c(0.2, 3, 20)
fraction_sort <- c("pico", "nano")

sample_type_filtered <- c(str_c(sample_type[1], fraction_filter, "um", sep = " "), 
    str_c(sample_type[2], fraction_filter, "um", sep = " "), str_c(sample_type[3], 
        fraction_sort, sep = " "))

# Export as specific data set as a phyloseq file

for (one_sample_type in sample_type[1:2]) {
    for (one_fraction in fraction_filter) {
        one_sample_type_filtered <- c(str_c(one_sample_type, one_fraction, "um", 
            sep = " "))
        ps[[one_sample_type_filtered]] <- ps[[one_sample_type]] %>% subset_samples(fraction_min == 
            one_fraction) %>% filter_taxa(function(x) sum(x) > 0, TRUE)
        
        print(glue("Phyloseq - {one_sample_type_filtered}"))
        print(ps[[one_sample_type_filtered]])
        cat("============================\n")
    }
    
}

Phyloseq - 18S filter 0.2 um
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 330 taxa and 42 samples ]
sample_data() Sample Data:       [ 42 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 330 taxa by 8 taxonomic ranks ]
============================
Phyloseq - 18S filter 3 um
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 374 taxa and 43 samples ]
sample_data() Sample Data:       [ 43 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 374 taxa by 8 taxonomic ranks ]
============================
Phyloseq - 18S filter 20 um
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 306 taxa and 35 samples ]
sample_data() Sample Data:       [ 35 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 306 taxa by 8 taxonomic ranks ]
============================
Phyloseq - 16S plastid 0.2 um
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 181 taxa and 28 samples ]
sample_data() Sample Data:       [ 28 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 181 taxa by 8 taxonomic ranks ]
============================
Phyloseq - 16S plastid 3 um
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 219 taxa and 41 samples ]
sample_data() Sample Data:       [ 41 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 219 taxa by 8 taxonomic ranks ]
============================
Phyloseq - 16S plastid 20 um
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 212 taxa and 31 samples ]
sample_data() Sample Data:       [ 31 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 212 taxa by 8 taxonomic ranks ]
============================

# Sorted samples

for (one_sample_type in sample_type[3]) {
    for (one_fraction in fraction_sort) {
        one_sample_type_filtered <- c(str_c(one_sample_type, one_fraction, sep = " "))
        ps[[one_sample_type_filtered]] <- ps[[one_sample_type]] %>% subset_samples(fraction_name == 
            one_fraction) %>% filter_taxa(function(x) sum(x) > 0, TRUE)
        
        print(glue("Phyloseq - {one_sample_type_filtered}"))
        print(ps[[one_sample_type_filtered]])
        cat("============================\n")
    }
    
}

Phyloseq - 18S sort pico
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 178 taxa and 27 samples ]
sample_data() Sample Data:       [ 27 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 178 taxa by 8 taxonomic ranks ]
============================
Phyloseq - 18S sort nano
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 124 taxa and 27 samples ]
sample_data() Sample Data:       [ 27 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 124 taxa by 8 taxonomic ranks ]
============================

sample_type_all <- c(sample_type, sample_type_filtered)

Normalize and transform to long form

for (one_sample_type in sample_type_all) {
    
    cat(one_sample_type)
    cat("\n============================")
    
    ps[[one_sample_type]] <- phyloseq_normalize_median(ps[[one_sample_type]])
    
    long[[one_sample_type]] <- phyloseq_transform_to_long(ps[[one_sample_type]])
    
    cat("\n")
    
}

18S filter
============================
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 701 taxa and 120 samples ]
sample_data() Sample Data:       [ 120 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 701 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  24221
16S plastid
============================
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 417 taxa and 100 samples ]
sample_data() Sample Data:       [ 100 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 417 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  33411
18S sort
============================
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 251 taxa and 54 samples ]
sample_data() Sample Data:       [ 54 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 251 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  25560
18S filter 0.2 um
============================
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 330 taxa and 42 samples ]
sample_data() Sample Data:       [ 42 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 330 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  17735
18S filter 3 um
============================
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 374 taxa and 43 samples ]
sample_data() Sample Data:       [ 43 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 374 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  24832
18S filter 20 um
============================
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 306 taxa and 35 samples ]
sample_data() Sample Data:       [ 35 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 306 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  33410
16S plastid 0.2 um
============================
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 181 taxa and 28 samples ]
sample_data() Sample Data:       [ 28 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 181 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  29616
16S plastid 3 um
============================
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 219 taxa and 41 samples ]
sample_data() Sample Data:       [ 41 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 219 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  35272
16S plastid 20 um
============================
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 212 taxa and 31 samples ]
sample_data() Sample Data:       [ 31 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 212 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  29330
18S sort pico
============================
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 178 taxa and 27 samples ]
sample_data() Sample Data:       [ 27 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 178 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  28592
18S sort nano
============================
========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 124 taxa and 27 samples ]
sample_data() Sample Data:       [ 27 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 124 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  19493

Treemaps

treemap_class <- list()
treemap_genus <- list()

for (one_sample_type in sample_type_all) {
    
    treemap_class[[one_sample_type]] <- phyloseq_long_treemap(long[[one_sample_type]], 
        division, class, str_c(one_sample_type, " - Class"))
    treemap_genus[[one_sample_type]] <- phyloseq_long_treemap(long[[one_sample_type]], 
        class, genus, str_c(one_sample_type, " - Genus"))
    
}

ASVs bargraphs

bargraph_asv <- list()
bargraph_species <- list()

for (one_sample_type in sample_type_all) {
    
    bargraph_asv[[one_sample_type]] <- phyloseq_long_bargraph(long[[one_sample_type]], 
        text_scaling = 0.75, n_bars = 20, use_asv = TRUE, title = one_sample_type)
    bargraph_species[[one_sample_type]] <- phyloseq_long_bargraph(long[[one_sample_type]], 
        text_scaling = 0.75, n_bars = 20, use_asv = FALSE, title = one_sample_type)
    
}

# Barplots per sample

bargraph_sample <- list()

for (one_sample_type in sample_type_all[4:11]) {
    
    bargraph_sample[[one_sample_type]] <- plot_bar(ps[[one_sample_type]], x = "sample_label", 
        fill = "class") + geom_bar(aes(color = class, fill = class), stat = "identity", 
        position = "stack") + ggtitle(str_c("Class level - ", one_sample_type)) + 
        theme(axis.text.y = element_text(size = 10)) + theme(axis.text.x = element_text(angle = 0, 
        hjust = 0.5)) + coord_flip() + scale_fill_viridis_d() + scale_color_viridis_d()
    print(bargraph_sample[[one_sample_type]])
}

Heatmaps (using the 10% most abundant taxa)

phyloseq_keep_abundant <- function(ps, contrib_min = 0.1) {
    total_per_sample <- max(sample_sums(ps))
    ps <- filter_taxa(ps, function(x) sum(x > total_per_sample * contrib_min) > 
        0, TRUE)
    ps <- phyloseq_normalize_median(ps)
    return(ps)
}

heatmap_class <- list()

for (one_sample_type in sample_type_all[4:11]) {
    
    ps_heat <- tax_glom(ps[[one_sample_type]], taxrank = "class") %>% phyloseq_keep_abundant(contrib_min = 0.1)
    
    heatmap_class[[one_sample_type]] <- plot_heatmap(ps_heat, method = "NMDS", 
        distance = "bray", taxa.label = "class", taxa.order = "division", sample.label = "sample_label", 
        sample.order = "sample_label", low = "beige", high = "red", na.value = "beige", 
        title = one_sample_type)
    
    # plot(heatmap(otu_table(ps_heat)))
    
    print(heatmap_class[[one_sample_type]])
    
}

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 7 taxa and 42 samples ]
sample_data() Sample Data:       [ 42 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 7 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  16751

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 7 taxa and 43 samples ]
sample_data() Sample Data:       [ 43 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 7 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  23530

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 5 taxa and 35 samples ]
sample_data() Sample Data:       [ 35 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 5 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  31889

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 6 taxa and 28 samples ]
sample_data() Sample Data:       [ 28 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 6 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  29445

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 5 taxa and 41 samples ]
sample_data() Sample Data:       [ 41 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 5 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  34645

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 4 taxa and 31 samples ]
sample_data() Sample Data:       [ 31 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 4 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  28471

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 6 taxa and 27 samples ]
sample_data() Sample Data:       [ 27 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 6 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  28038

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 5 taxa and 27 samples ]
sample_data() Sample Data:       [ 27 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 5 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  19227

heatmap_genus <- list()

for (one_sample_type in sample_type_all[4:11]) {
    
    ps_heat <- tax_glom(ps[[one_sample_type]], taxrank = "genus") %>% phyloseq_keep_abundant(contrib_min = 0.1)
    
    heatmap_genus[[one_sample_type]] <- plot_heatmap(ps_heat, method = "NMDS", 
        distance = "bray", taxa.label = "genus", taxa.order = "division", sample.label = "sample_label", 
        sample.order = "sample_label", low = "beige", high = "red", na.value = "beige", 
        title = one_sample_type)
    
    print(heatmap_genus[[one_sample_type]])
    
}

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 13 taxa and 42 samples ]
sample_data() Sample Data:       [ 42 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 13 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  15651

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 6 taxa and 43 samples ]
sample_data() Sample Data:       [ 43 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 6 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  20393

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 10 taxa and 35 samples ]
sample_data() Sample Data:       [ 35 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 10 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  30906

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 6 taxa and 28 samples ]
sample_data() Sample Data:       [ 28 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 6 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  28612

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 7 taxa and 41 samples ]
sample_data() Sample Data:       [ 41 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 7 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  33791

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 5 taxa and 31 samples ]
sample_data() Sample Data:       [ 31 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 5 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  27659

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 11 taxa and 27 samples ]
sample_data() Sample Data:       [ 27 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 11 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  27435

========== 
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 6 taxa and 27 samples ]
sample_data() Sample Data:       [ 27 samples by 61 sample variables ]
tax_table()   Taxonomy Table:    [ 6 taxa by 8 taxonomic ranks ]

==========
The median number of reads used for normalization is  17950

Figures

treemap_18S_filt <- cowplot::plot_grid(treemap_class[["18S filter 0.2 um"]]$gg, 
    treemap_class[["18S filter 3 um"]]$gg, treemap_class[["18S filter 20 um"]]$gg, 
    ncol = 1)

treemap_18S_sort <- cowplot::plot_grid(treemap_class[["18S sort pico"]]$gg, 
    treemap_class[["18S sort nano"]]$gg, ncol = 1)

treemap_16S_filt <- cowplot::plot_grid(treemap_class[["16S plastid 0.2 um"]]$gg, 
    treemap_class[["16S plastid 3 um"]]$gg, treemap_class[["16S plastid 20 um"]]$gg, 
    ncol = 1)

fig_2 <- cowplot::plot_grid(bargraph_species[["18S filter"]]$gg, treemap_18S_filt, 
    bargraph_species[["18S sort"]]$gg, treemap_18S_sort, bargraph_species[["16S plastid"]]$gg, 
    treemap_16S_filt, ncol = 2, rel_widths = c(2, 1))

fig_2